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Abstract 

We derive general kinetic and hydrodynamic models of chemotactic aggregation that 
describe certain features of the morphogenesis of biological colonies (like bacteria, amoe- 
bae, endothelial cells or social insects). Starting from a stochastic model defined in terms 
of N coupled Langevin equations, we derive a nonlinear mean field Fokker-Planck equa- 
tion governing the evolution of the distribution function of the system in phase space. By 
taking the successive moments of this kinetic equation and using a local thermodynamic 
equilibrium condition, we derive a set of hydrodynamic equations involving a damping 
term. In the limit of small frictions, we obtain a hyperbolic model describing the for- 
mation of network patterns (filaments) and in the limit of strong frictions we obtain a 
parabolic model which is a generalization of the standard Keller-Segel model describing 
the formation of clusters (clumps). Our approach connects and generalizes several mod- 
els introduced in the chemotactic literature. We discuss the analogy between bacterial 
colonies and self-gravitating systems and between the chemotactic collapse and the gravi- 
tational collapse (Jeans instability). We also show that the basic equations of chemotaxis 
are similar to nonlinear mean field Fokker-Planck equations so that a notion of effective 
generalized thermodynamics can be developed. 



Key words: Nonlinear mean field Fokker-Planck equations, generalized thermodynam- 
ics, chemotaxis, gravity, long-range interactions 
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1 Introduction 



In many fields of physical sciences, one is confronted with the description of the evolution of a 
system of particles which self-consistently attract each other over large distances [U This 
is the case in biology, for example, in relation with the process of chemotaxis [3]. Chemotaxis 
explains the spontaneous self-organization of biological cells (bacteria, amoebae, endothelial 
cells,...) or even insects (like ants) due to the long-range attraction of a chemical (pheromone, 
smell, food,...) produced by the organisms themselves. The chemotactic aggregation of biolog- 
ical populations is usually studied in terms of the Keller-Segel model consisting in two coupled 
differential equations [I]: 

^ = jD ,Ap- x V-(pVc), (1) 

dc 

— = -k( c )c + h(c)p + D c Ac, (2) 

which describe the evolution of the concentration p(r, t) of the biological organisms and of the 
secreted chemical c(r, t). The Keller-Segel model is a parabolic model where the evolution of 
the concentration of the biological organisms is governed by a drift-diffusion equation ([1]). The 
diffusion models the erratic motion of the particles (like in Brownian theory) and the drift 
term models a systematic motion along the gradient of concentration of the secreted chemical. 
When x > the cells are attracted in regions of high concentration while for \ < they are 
repelled from the regions of high concentration (in that case the chemical acts as a poison). 
The evolution of the secreted chemical is described by a diffusion equation (J2J) involving terms 
of source and degradation: the chemical is produced by the organisms at a rate h(c) and 
it is degraded at a rate fc(c). For x > 0, the Keller-Segel model is able to reproduce the 
chemotactic aggregation (collapse) of biological populations when the attractive drift term 
Xp||Vc|| overcomes the diffusive term D^||Vp|| above a critical mass M c [5]. This is similar to 
the gravitational collapse of self-gravitating Brownian particles, described by the Smoluchowski- 
Poisson system, below a critical temperature T c [6] (see [H [8] for a detailed discussion of 
the analogy between the Keller-Segel model and the Smoluchowski-Poisson system). These 
parabolic models ultimately lead to the formation of Dirac peaks [9j [10] . 

Some regularizations of the Keller-Segel model have been introduced. They have the form 
of generalized drift-diffusion equations [TTJ [7] : 

§ = V-( X (Vp-pVc)), (3) 

dc 

— = -k( c )c + h(c)p + D c Ac, (4) 

in which the simple diffusion term -D*Vp in Eq. (CEJ) is replaced by a more general "pressure 
term" Vp(p). The drift-diffusion equation ([3]) is similar to a generalized Smoluchowski equation 
[T2| [T3] 3- By adapting the barotropic equation of state p(p), one can obtain regularized 
chemotactic models preventing the density from reaching infinitely large values and forming 
singularities [14] . In that case, the Dirac peaks (clumps) are replaced by smoother density 
profiles (aggregates). The dynamical evolution of the regularized model ©-(BJ generically 
leads to the formation of A/"(£) round aggregates which progressively merge until only one big 
aggregate remains at the end. 



1 The possibility that the diffusion coefficient D* and the chemotactic sensitivity x depend on the density of 
cells p and on the density of chemical c is considered in the primitive model of Keller & Segel [4] (see Appendix 
IE.4[) . The very much studied model (HI)-© is a simplification of the primitive Keller-Segel model where the 
coefficients D* and x are assumed constant. 
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However, recent experiments of in vitro formation of blood vessels show that cells randomly 
spread on a gel matrix autonomously organize to form a connected vascular network that is 
interpreted as the beginning of a vasculature [15]. This phenomenon is responsible of angio- 
genesis, a major actor for the growth of tumors. These networks cannot be explained by the 
parabolic models flH)-(|ID that l ea, d to pointwise blow-up or round aggregates. However, they 
can be recovered by hyperbolic models that lead to the formation of networks patterns that are 
in good agreement with experimental results. These models take into account inertial effects 
and they have the form of hydrodynamic equations [151 : 



^ + V-(pu) = 0, (5) 

-£ + (u-V)u=— Vp + Vc, (6) 
ot p 

Oc 

— = -k( c )c + h(c)p + D c Ac (7) 
at 

The inertial term models cells directional persistence and the general density dependent pressure 
term — Vp(p) can take into account the fact that the cells do not interpenetrate. In these models, 
the particles concentrate on lines or filaments. These structures share some analogies with the 
formation of ants' networks (due to the attraction of a pheromonal substance) and with the 
large-scale structures in the universe that are described by similar hydrodynamic (hyperbolic) 
equations: the so-called Euler-Poisson system. The similarities between the networks observed 
in astrophysics (see Figs 10-11 of [16]) and biology (see Figs 1-2 of [15] ) are striking. 

The above-mentioned parabolic and hyperbolic models are continuous models which describe 
the evolution of a smooth density field p(r, t) and, in the case of hyperbolic models, a smooth 
velocity field u(r, t). In this paper, we propose a kinetic derivation of these models starting from 
a microscopic description of the dynamics of the biological population. We introduce stochastic 
equations for the motion of each individual and, implementing a mean field approximation, we 
obtain the corresponding Fokker-Planck equation governing the evolution of the distribution 
function f(r,v,t) of the system in phase space. The stochastic Langevin equations involve 
a friction force, an effective force due to the chemotactic attraction of the chemical and a 
random force (noise) whose strength can depend on the local distribution function itself. This 
dependence can take into account microscopic constraints ("hidden constraints") that affect 
the dynamics of the particles at small scales. We can have (i) close packing effects (like finite 
size effects, excluded volume constraints, steric hindrance...) that forbid the interpenetration 
of the particles and prevent the system from reaching arbitrarily high densities [2] and (ii) 
nonextensivity effects that alter the usual random walk and lead to anomalous diffusion and 
non-ergodic behaviour [T71[T8]. The resulting generalized stochastic equations lead to nonlinear 
mean field Fokker-Planck equations similar to those occurring in the context of generalized 
thermodynamics [191 E2l HI]- Therefore, as first noticed in [11], the chemotaxis of biological 
populations can be a physical system where a notion of effective generalized thermodynamics 
applies. By taking the successive moments of these generalized Fokker-Planck equations and 
using a local thermodynamic equilibrium condition, we derive a closed set of hydrodynamic 
equations 

^ + V-(pu) = 0, (8) 

Ji + (u-V)u = — Vp + Vc-Cpu, (9) 
ot p 

dc 

= -k( c )c + h(c)p + D c Ac, (10) 

at 
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involving a friction term — £u [H] ED]- For £ = 0, we recover the hyperbolic model (JS])-(J7J) pro- 
posed by Gamba et al. [15] to model vasculogenesis and for £ — > +00 we obtain the parabolic 
model (jSI)-© generalizing the Keller-Segel model ([I])-®. We discuss the analogy between bac- 
terial colonies and self-gravitating systems and between the chemotactic collapse and the "grav- 
itational collapse" (Jeans instability) [TJ [HJ [2D] • Indeed, the kinetic and hydrodynamic models 
of biological populations derived in this paper are similar to those describing self-gravitating 
systems [HI [21]. Therefore, our approach connects various topics studied by different commu- 
nities: systems with long-range interactions pQ, nonlinear mean field Fokker-Planck equations 
and generalized thermodynamics [19J, self-gravitating systems [22] and chemotaxis [I]. 

2 A stochastic model of chemotactic aggregation 
2.1 Generalized Langevin equations 

We shall introduce a model of chemotactic aggregation generalizing the Keller-Segel model (pQ)- 
([2]). For biological systems, the number of constituents is not necessarily large so that it may 
be relevant to return to a "corpuscular" description of the system and introduce an equation 
of motion for each particle (= cell). This type of "microscopic" approach has been previously 
considered by Schweitzer & Schimansky-Geier [23] , Stevens [23] and Newman & Grima [25] (see 
also related work by the authors 0). They describe the motion of individuals by a stochastic 
equation of the form 

^ = xV a c+y2aR a (t). (11) 

The first term in the r.h.s. is the chemotactic drift to which the particles are submitted (the 
coefficient x plays the role of a mobility). The second term is a stochastic term where R Q (i) is a 
white noise satisfying (R a (t)) = and (R ita (t)Rj t p(t')) = 8ij5 a> p5{t—t') (where a = 1, N refer 
to the particles and i = 1, d to the space coordinates) and is a diffusion coefficient. The 
diffusion, that is observed for several biological organisms, can have different origins depending 
on the system under consideration [27]. In the case of small organisms moving in a fluid 
(matrigel), it can be due to the repeated impact of the molecules of the fluid on the particles 
like in ordinary Brownian motion for colloidal suspensions. In other cases, it can be due to the 
properties of motion of the particles themselves. For example, bacteria like Escherichia coli 
are equipped with flagella and are self-propelled. When rotated counterclockwise, the flagella 
act as a propellor and the bacterium moves along straight line ("run"). Suddenly, the flagella 
rotate clockwise and the bacterium stops to choose a new direction at random ("tumble"). It 
continues in that direction for a while until the next tumble. Therefore, the bacteria experience 
a random motion of their own. At the simplest level of description, this motion can be modelled 
by a stochastic term like in Eq. (II ip . These stochastic equations describe the motion of each 
of the N particles of the colony. As indicated above, N is not necessarily large so it may 

2 Chavanis et al. [7J [H] studied stochastic models of Brownian particles with long-range interactions 
where the motion of the individuals is described by stochastic equations of the form coupled by a binary 
potential of interaction u{\ti — |) instead of the more complicated potential c, solution of Eq. (Tl2|) . depending 
on the past history of the system (memory terms are, however, considered in [71ll3j). These models describe, for 
example, self-gravitating Brownian particles |21| and a simplified chemotactic model where Eq. (|12p is replaced 
by Ac— fcpC = — A^^ =1 S(r — r a (t)). This simplification is valid in a limit of large diffusivity of the chemical (see 
Appendix [C| so that dc/dt can be neglected. In the mean field approximation valid for N — ► +00 in a proper 
thermodynamic limit :26J, these models reduce to the Smoluchowski-Poisson system [6] or to the Keller-Segel 
model ((T|)- ^ of chemotaxis where Eq. @ is replaced by Ac — fcgc = — Xp [14] , Note that microscopic models 
yielding the regularized Keller-Segel model (J3])— (J4j) have also been introduced in [14] . 
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be of interest to treat the bacterial colony as a discrete system of particles. By contrast, the 
chemical that is secreted is usually described as a continuous field. Therefore, the evolution of 
the concentration of the chemical is governed by an equation of the form 



dc 



dt 



-kc + D c Ac + hJ2 5 ( r -r a (t)). (12) 



a=l 



Equations ( JTTl) - (fl2l) have been studied in [2H EH [25J [26]. In the mean-field approximation, 
they return the usual Keller-Segel model (JT])-(j2J). 

We shall generalize the model (111I) - (TI21) in two respects. First of all, there exists biological 
systems for which the inertia of the particles has to be taken into account [15J. This "inertia" 
means that they do not respond immediately to the chemotactic drift. We propose therefore 
to describe the motion of each individual of the biological population by a stochastic equation 
of the form 

(/V " -^v a + V a c + V2DR a (t). (13) 



dt dt 

The first term in the r.h.s. is an "effective" friction force, the second term is a force that 
models the chemotactic attraction due to the chemical c and the last term is a random force. 
The friction force takes into account the fact that the velocity of the particles has the tendency 
to be directed along the concentration gradient Vc. This is exactly the case when £ — > +oo. In 
this strong friction limit, we recover the overdamped model (11 II) with x — l/£ an d -D* = -D/£ 2 - 
For finite values of £, the velocity will take a (relaxation) time r ~ to get aligned with 
the concentration gradient. This is how an inertial effect is introduced in the model. The 
term — £v can also represent a physical friction of the organisms against a fixed matrigel. 
The stochastic equations (fl3"j) . coupled self-consistently to the field equation ffl2|) . describe the 
motion of each of the N particles of the colony. This completely discrete model of chemotactic 
aggregation, which takes into account statistical correlations between the particles, is developed 
in Appendix |A] An exact equation for the single-cell probability distribution is derived and 
it is shown precisely how a mean field approximation can be implemented in the theory. In 
the mean-field approximation, passing to a hydrodynamical description, we obtain the model 
(1H|)- (|T0|) with a linear (isothermal) equation of state p = pT e ff — Dp/£. 

In order to describe more general situations where the equation of state p = p(p) is nonlinear, 
we shall consider a generalized class of stochastic equations where the diffusion coefficient 
explicitly depends on the distribution function /(r, v,t) of particles in phase space. These 
generalized stochastic equations have been introduced in [281 El US H] . They are associated 
with nonlinear Fokker-Planck equations and lead to a notion of "generalized thermodynamics" . 
The possibility to apply this type of equations to the chemotactic problem was proposed in 
[TTj . In order to simplify the formalism, we shall make a mean-field approximation since the 
start and describe the motion of each individual of the biological population by a stochastic 
equation of the form 

^ = + V S(c) + v^Rfi), (14) 
where the mean field force is determined by the equation 

dc 

— = -k( c )c + h(c)p + D c Ac, (15) 

where p(r, t) = (^^(r — r^i))) is the smooth local density of cells (the brackets denote an 
average over the noise). This amounts to replacing the exact density p ex (j, t) = 5(r — rj(t)) 
in Eq. ( TT21) by the smooth density p(r,t) = (X]j^( r — r i(^)))- This is how the mean-field 



5 



approximation is introduced in the model. For sake of generality, we have allowed the coefficients 
k and h to depend on the concentration c and we have written the chemotactic force as the 
gradient of a function S(c) of the concentration. Equations (jT4"|) - (IT5|) will be our starting point 
for the chemotactic problem. This model, or the discrete model (1131) . is similar to the model of 
self-gravitating Brownian particles introduced in [SIEZIEI]- The main difference, beyond the 
context, is that the Poisson equation in gravity is replaced by the more general field equation 

(USD. 



2.2 Nonlinear mean field Fokker-Planck equations 

We shall now derive kinetic and hydrodynamic equations associated with the stochastic model 
f fT4|) - f|T5|) . Using standard methods of Brownian theory, we find that the evolution of the 
distribution function f(r, v,t) of the system is described by a generalized mean field Fokker- 
Planck equation of the form 



^ + v .^ + V 5(c)-^ = — 
dt dr dv dv 



d_ 
dv 



(D(f)f) + £fv 



(16) 



We have considered a relatively general class of stochastic processes ( !T4l) where the diffusion 
coefficient D(f) can depend on the distribution function. This dependence can take into account 
microscopic constraints of various origin ("hidden" constraints) that affect the dynamics at 
small scales. These generalized stochastic processes are associated with a notion of generalized 
thermodynamics [19] (see also Appendix [E]). The case where D(f) is a power law has been first 
considered by Borland [28] in connection with Tsallis generalized thermodynamics [29]. The 
case where D(f) is arbitrary has been considered by Chavanis [11J in connection with nonlinear 
Fokker-Planck (NFP) equations associated with more general forms of entropic functionals than 
the Tsallis entropy. If we write the diffusion coefficient in the form D(f) = Df[C(f)/f]', where 
C(f) is a convex function (i.e. C" > 0), the nonlinear Fokker-Planck equation (Tl6l) can be 
rewritten [TT] : 

dj_ 8£ 8£ _ d_ 

dt dr dv dv 

This type of nonlinear Fokker-Planck equations can also be derived from a master equation 
and a Kramers-Moyal expansion by allowing the transition probabilities to depend on the 
occupation number in the initial and arrival states as done in Kaniadakis [30]. The case of 
a normal diffusion D(f) = D corresponds to C(f) = /In/ leading to the usual Kramers 
equation fl66l) associated with the Boltzmann statistics. The choice D(f) = Df q ~ x corresponds 
to C(f) = (f q — f)/(q— 1) leading to the polytropic Kramers equation [2H dH [T7J Q2] associated 
with the Tsallis statistics (these examples will be worked out explicitly in Sec. 12. 4p . 

The stationary solution of Eq. (TTT1) . which cancels both the "collision" term (r.h.s.) and 
the advection term (l.h.s.), is given by (see Appendix [El): 



DfC"(f)^ + £fv 



(17) 



C'(f, 



eq) 



V 

~2 



S(c) 



— oi, 



where (3 = 1/T e fj is an effective inverse temperature satisfying a generalized Einstein relation 
j3 — £/D and a is a constant of integration. Since C is convex, the above relation can be 
reversed. Then, we find that / eg (r, v) = F(f3e + a) = / eg (e) where F(x) = (C") _1 (— x). Since 

and (3 > 0, we find that f eq (e) is a decreasing function of the individual 



-P/c'U 

,2 



eq ) 



energy e = v 2 /2 — S(c). The case of normal diffusion D(f) = D leads to the Maxwell-Boltzmann 
distribution f eq = Ae~ l3t . The case of anomalous diffusion D(f) = Df^ 1 leads to the Tsallis 



distribution f eq = [A — (3{q — l)e/q] 



1/(9-1) 
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2.3 Damped hydrodynamic equations 

We shall now derive the moments equations issued from the generalized Fokker-Planck equation 
( TT7|) . Defining the density and the local velocity by 

p = J fdv, pu= J fvdv, (19) 

and integrating Eq. (ITTj) on velocity, we get the continuity equation 

^ + V • (pu) = 0. (20) 

Next, multiplying Eq. (I17p by v and integrating on velocity, we obtain 

d d dP~ ' dc 

m (pm) + Q^Ows) = + pS'(c)^ - {pa*, (21) 

where we have defined the "pressure" tensor 

Pij = [ fWiWj dv, (22) 



where w = v — u is the relative velocity. Using the continuity equation, Eq. (JzTj) can be 
rewritten 

By taking the successive moments of the velocity, we can obtain a hierarchy of hydrody- 
namic equations. Each equation of the hierarchy involves the moment of next order. If we 
are sufficiently close to equilibrium, it makes sense to close the hierarchy of equations by using 
a condition of local thermodynamic equilibrium (L.T.E.). We shall thus evaluate the pres- 
sure tensor Eq. fl22|) with the distribution function fL.T.E( r , v, t) defined by the relation (see 
Appendix EI]): 



C'Ul.t.e) = -P 



w 2 

Y + A(r,t) 



(24) 



The function A(r, t) is implicitly related to the density by writing 

p(r,t) = J f L .T.E.dv = p[\(r,t)}. (25) 

Using the condition Eq. (1241) of local thermodynamic equilibrium, the pressure tensor Eq. ( 1221) 
can be written = pSij with 

p(r,t) = - d J f L . T . E w 2 dw=p[\(v,t)}. (26) 

The pressure is a function p = p(p) of the density which is entirely specified by the function 
C(f), by eliminating A from the relations Eq. (1251) and Eq. (|26|) . This defines a barotropic 
gas. For example, the case of normal diffusion D(f) = D leads to a Maxwellian distribution 
fh.T.E. = ((3/2Tc) d / 2 p(r,t)exp(—(3w 2 /2) and a linear equation of state p = pT e ff as for an 
isothermal gas. The case of anomalous diffusion D(f) = Df 1 ^ 1 leads to a Tsallis distribution 
Jl.t.e. = [p( r )^) — /3((q — ]-)/(l)w 2 /2] 1 ^ q ^ and a power law equation of state p = Kp 1 (with 
7 = 1 + \jn and n = d/2 + l/(q — 1)) as for a polytropic gas [T7]. Substituting the result P^ = 
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p(p)Sij in Eq. (j23J) and collecting the other constitutive equations, we obtain a hydro dynamic 
model of the form: 

^ + V • (pu) = 0, (27) 

+ (u ■ V)u = — Vp + V5(c) - £u, (28) 
or p 

<9c 

_ = -k(c)c + h(c)p + D c Ac. (29) 
at 

The damped barotropic Euler equations (l2Tj) - (l29l are interesting as they connect hyperbolic 
models to parabolic models [11]. For £ = 0, we recover the hydrodynamic model (E])- (JTJ) 
introduced by Gamba et al. [15j IE Alternatively, for £ — > +oo, we can formally neglect the 
inertial term in Eq. (1281 so that the velocity field is given by 

u = -^(Vp - pS'(c)Vc) + 0(r 2 )- (30) 

Substituting this drift term in the continuity equation (1271) . we obtain the drift- diffusion equa- 
tion 



at 



~(Vp-pS'(c)Vc) 



(31) 



which is a generalization of the Keller-Segel model. The usual Keller-Segel model is recovered 
for a normal diffusion D(f) = D leading to p = pT e ff = Dp/£. We make the link with Eq. 
(PQ) by setting = D/£ 2 , \ — l/£ an d ^(c) = c - The case of a polytropic equation of state 
p = i\p 7 associated to the Tsallis statistics has been studied in [T71 [18] . 

It should be stressed that the damped Euler equations (I2"71 - (j2"9l remain heuristic because 
their derivation is based on the Local Thermodynamic Equilibrium (L.T.E.) condition (I24j) 
which is not rigorously justified. However, using a Chapman- Enskog expansion, it is shown 
in [33] that the generalized Smoluchowski equation (I3T1) is exact in the limit £ — > +oo (or, 
equivalently, for times t ^> £ _1 ). The generalized Smoluchowski equation can also be obtained 
from the moments equations of the generalized Kramers equation by closing the hierarchy in 
the limit £ — > +oo (see [13] and Appendix IB"]) . 



2.4 A generalized Keller-Segel model 

In this section, we discuss an explicit example to illustrate our general formalism. We consider 
a stochastic process of the form 

^ = _£v + V5(c) + </2DfW*R(t), (32) 

where the mean field force is determined by Eq. (fT5|) . Comparing with Eq. (fl4j) . we find that 
the diffusion coefficient is given by D(f) = Df q ~ x . As indicated previously, this leads to a 
situation of anomalous diffusion related to the Tsallis statistics [28]. For q = 1, we recover 

3 In fact, the derivation of the damped barotropic equations pT)) -(J!?!) f from the kinetic theory developed in 
this section implicitly assumes that the friction force — £u is sufficiently large. Therefore, the limit £ — * is not 
really justified. Indeed, when D, £ — > 0, Eq. (fT6|) reduces to the Vlasov equation and the L.T.E condition (|24|) is 
not realized or takes a long time to establish itself. In that collisionless regime, the system can undergo a process 
of "violent relaxation" driven by mean-field effects, like in astrophysics for the Vlasov- Poisson system [32]. In 
Appendix [Dj we present an alternative kinetic theory, based on other assumptions, leading to the hyperbolic 
model ©-© proposed by Gamba et al. [T3] , 
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the standard Brownian model with a constant diffusion coefficient, corresponding to a pure 
random walk (see Appendix |A|) . In that case, the sizes of the random kicks are uniform and 
do not depend on where the particle happens to be. For q ^ 1, the size of the random kicks 
changes, depending on the phase-space distribution of the particles around the "test" particle. 
A particle which is in a state (r, v) that is highly populated [large f(r, v, t)} will tend to have 
larger kicks if q > 1 and smaller kicks if q < 1. Since the microscopies depends on the actual 
density in phase space, this creates a bias in the ergodic behavior of the system. The nonlinear 
Fokker-Planck equation associated with the stochastic process (1321) is the polytropic Kramers 
equation 

9 / +v .^ + VS (c).^^.f^ +f /vY (33) 



dt dr <9v <9v \ <9v 

For q — 1, we recover the classical Kramers equation (1661) which can be deduced from the 
iV-body stochastic process (fT3l) by using a mean field approximation (see Appendix |A|). For 
q 7^ 1, the NFP equation (133]) can be deduced from the mean field stochastic process (l32l) . The 
associated Lyapunov functional (11001) is explicitly given by 

F[f] = J fjdrdv + i- J [D c (Vcf + kc 2 ] dr- J pedr + I^L J (f q - f)drdv, (34) 

where T e jf = 1/(3 = D /£. It can be viewed as an effective generalized free energy F = E—T e ffS 
associated with the Tsallis entropy S q = —(l/(q — 1)) J (f q — f)drdw. It satisfies F < which 
is the version of the if -theorem in the canonical ensemble where the temperature T e ff is fixed 
instead of the energy. The steady state of the NFP equation (!33j) is the polytropic (Tsallis) 
distribution 



eq 



H e , (35) 

1 

where e = v 2 /2 — S(c) is the energy per particle. This distribution extremizes the free energy 
( T34|) at fixed mass. Furthermore, it is linearly dynamically stable with respect to the NFP 
equation (1331) if, and only if, it is a minimum of F at fixed M [11J. The index n of the 
polytrope (defined by Eq. (f4"T|) below) is related to the parameter q by the relation 

d 1 

n = - + — . (36) 

The isothermal (Boltzmann) distribution function f eq = Ae~^ e , corresponding to normal diffu- 
sion, is recovered in the limit q — > 1, i.e. n — > +oo. In the following, we shall consider q > 
so that C is convex. Since (3 > 0, this implies that /(e) is decreasing (see Sec. 12. 2p . There are 
two cases to consider. For q > 1, i.e. n > d/2, the distribution (1351) can be written 

/ = A(e m - e) 1 ^" 1 ) (easel) (37) 

where A = \(3{q — \)/q} 1 ^ q ~ 1 ' > and e m = qfi/[P(q — 1)]. It has a compact support since / is 
defined only for e < e m . For e > e m , we set / = (for q = +oo, i.e. n = d/2, f is the Heaviside 
function). For q < 1 the distribution can be written 

/ = A(e + e)- 1/{1 - q) (case 2) (38) 

where A = [(3(1 — q) / g] -1 /* 1 -?) and eo = qn/[(3(l — q)\. It is defined for all energies. For large 
velocities, it behaves like / ~ v 2n ~ d _ Therefore, the density and the pressure are finite only for 
n < — 1, i.e. d/(d + 2) < q < 1. Therefore the range of allowed parameters are 

q > 1, n > — (case 1) (39) 



9 



d 



d + 2 



<q<l, 



n < — 1 (case 2) 



(40) 



From Eq. ( 1351) . or the alternative forms ( 1371) and ( 1381) . we can compute the density p = f fdv = 
p[c(r)] and the pressure p = (1/d) f fv 2 dw = p[c(r)] at equilibrium (see Sec. 3.6 of [34J for 
details). Then, eliminating the concentration c(r) between these expressions, we obtain the 
polytropic equation of state 

1 



p = Kp" 1 , 7 = 1 + 
For n > d/2 (case 1), the polytropic constant is 



n 



K 



n+ 1 



AS d 2? 



d_ 1 T(d/2)T(l-d/2 + n) 



and for n < — 1 (case 2), we have 



K 



n + 1 



AS d 2* 



T(l + n) 



-ir(d/2)r- 



-l/n 



T(d/2-n) 



-l/n 



(41) 



(42) 



(43) 



where T(x) is the Gamma function. For q — 1, the polytropic equation of state (|4T1) reduces 
to the "isothermal" (linear) equation of state p = pT e ff. These equations of state have been 
obtained from the distribution ( 1351) valid at equilibrium. However, the same equations of 
state are obtained from the L.T.E condition ( 1241) . using Eqs. ( 1251) and ( 1261) and C(/) = 
(f q — f) I (q — 1), or from the distribution function (176]) valid in the strong friction limit £ — > +oo, 
using Eqs. (1771) and (1781) . Then, the damped hydrodynamic equations corresponding to the 
model (132]) are Eqs. (!271) - (l29l) with the equation of state (|4TI) . On the other hand, in the strong 
friction limit, the drift-diffusion equation (13T1) reduces to the polytropic Smoluchowski equation 



dp 
di 



GFTVp 7 -pS'(c)Vc) 



The associated Lyapunov functional (11271) is explicitly given by 



F[P\ 



(44) 



(45) 



It can be viewed as an effective generalized free energy F = E — KS associated with the Tsallis 
entropy S 1 = —(1/(7 — 1)) J(p 7 — p)dv where K plays the role of an effective temperature 
and 7 the role of the Tsallis g-parameter. It satisfies F < which is the version of the H- 
theorem in the canonical ensemble where the "polytropic temperature" K is fixed instead of 
the energy. The steady state of the generalized Smoluchowski equation (jUj) is the polytropic 
(Tsallis) distribution in physical space 



A + 



7-1 



S(c) 



7-1 



(46) 



This distribution extremizes the free energy (1451) at fixed mass. Furthermore, it is linearly 
dynamically stable with respect to the generalized Smoluchowski equation (jUJ) if and only if it 
is a minimum of F at fixed M. Of course, Eq. (j46H can be deduced from Eq. (135|) . Equation (|44l) 
can be interpreted as a NFP equation of the form considered in [11] with C(p) = (p 7 — p)/(j — 1)- 
It can be obtained directly from the generalized stochastic process (see Appendix |E.4p 



dr 

dt 



P VS{c) 




(47) 
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This stochastic process generalizes Eq. ffTTl) and provides another justification of the polytropic 
Smoluchowski equation (jHJ). The interpretation of the noise in Eq. (14T1) . depending on the 
density p(r, t), is similar to that given after Eq. (1321) although the process takes place in physical 
space instead of phase space. 

Finally, we note the remarkable feature that a polytropic distribution f eq (e) [see Eq. (135"]) ] 
with index q in phase space yields a polytropic distribution p eq (c) [see Eq. ( 115]) ] with index 
7 in physical space (using Eqs. fl36l) and (j4T|) . the indices are related to each other by 7 = 
1 + 2(q — l)/[2 + d(q — 1)]). In this sense, polytropic (Tsallis) distributions are "stable" laws. 
Apparently, these are the only ones enjoying the property that f e q{t) and p eq {c) have the same 
form. By comparing Eqs. ( 1351) and ( 146]) or Eqs. ( 1341) and ( )45j) . we note that K plays the same 
role in physical space as the temperature T e ff = 1/(3 in phase space. It is sometimes called a 
"polytropic temperature". We also note that for q > 1, we have 7 > 1 so that the model is 
hyper-diffusive in phase space and physical space. For q < 1, we have 7 < 1 so that the model 
is sub-diffusive in phase space and physical space. For q = 1, we have 7 = 1 so that Eq. (|44l) 
reduces to the ordinary Smoluchowski equation. This yields the standard Keller-Segel model 
[I]. For q ^ 1, or 7 ^ 1, we obtain a generalized Keller-Segel model that can take into account 
non-ideal effects giving rise to anomalous diffusion. This model has been studied in [T7J \TE\ 155]. 



2.5 Numerical simulations 

In this section, we restrict ourselves to normal diffusion and we present numerical simulations 
of the iV-body system f[T5]) - f[T2]) in a simplified setting. In the equation f[T2]) for the evolution 
of the secreted chemical, we set h = XD C and consider a limit of large diffusivity D c — > +00 
with A ~ 1. In that limit, the temporal derivative dc/dt and the degradation term —kc can be 
neglected (see Appendix |C]) . The model (fT3l) - (fr2|) becomes 

" " -^ a + V a c + V2DR a (t), (48) 



dt 

Ac= -X(p ex -p), (49) 

where p ex (r,t) = ^^=1 — r a(^)) is the exact density and p = N/V is the average density 
over the entire domain. When D = £ = 0, we obtain 

dv, 



dt V « C ' (50) 

Ac= -X(p ex -p). (51) 

These equations are similar to the Newton equations for a self-gravitating system where — c 
plays the role of the gravitational potential $ and A the role of the gravitational constant SjG 
(Sd is the surface of a unit sphere in ^-dimensions) . In cosmology, when we take into account the 
expansion of the universe and work in a comoving frame, the usual Poisson equation A$ = 4ttGp 
is replaced by an equation of the form A0 = 47rGa(t) 2 [p(x, t) — Pbif)] where the density p(x, t) 
is replaced by the deviation p(x, t) — p b (t) to the mean density [SB]. Furthermore, the equations 
of motion read (d/dt)(ma 2 x) = — mV0. If we consider timescales over which the variation 
of the scale factor a{t) can be neglected, the equations of motion become isomorphic to Eqs. 
(150|) - (!5T|) . Therefore, there exists interesting analogies between the process of chemotaxis in 
biology and the dynamics of self-gravitating systems. 

However, in biology, we are rather in a limit where inertial effects are negligible or weak so 
that the (effective) friction coefficient £ is relatively large. Therefore, the process of chemotaxis 
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Figure 1: iV-body simulation of the inertial model fl50l) - fl5Tl) for iV = 32768 particles in a periodic 
box. The dynamical equations are isomorphic to the Newton equations for a self-gravitating 
system in cosmology They show the formation of a network pattern with a filamentary struc- 
ture. In biology, this corresponds to the beginning of a vasculature. 




Figure 2: iV-body simulation of the overdamped model (!52l - (l53l) for iV = 32768 particles in 
a periodic box. The dynamical equations are isomorphic to the Langevin equations for a self- 
gravitating Brownian gas in an overdamped limit. They lead to point-wise blow-up forming 
ultimately Dirac peaks. In biology, this corresponds to a chemotactic collapse. 
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in biology is more directly analogous to the dynamics of self-gravitating Brownian particles 
studied in [6]. If we consider the strong friction limit £ — > +00, we obtain 

v a = ^V a c + y/2D^R a (t), (52) 

Ac = -\{p ex - p). (53) 

These equations are similar to those describing a self-gravitating Brownian gas in an over- 
damped limit [2Tj . 

The inertial model (]50]) - (|5"T]) leads to fluid models of a form related to Eqs. ©-© that are 
hyperbolic. Alternatively, the overdamped model fl52|) - fl53|) leads to drift-diffusion equations 
of the form (JT|)- (J2J) that are parabolic. As discussed in the Introduction, parabolic models 
are known to lead to pointwise blow-up while hyperbolic models generate network patterns 
[T5] 137] . We have performed direct numerical simulations of the iV-body systems ( l50i) -( l5Ti) 
and (152]) - (153]) in a periodic domain starting from a statistically homogeneous distribution of 
particles. In order to simulate a large number of particles, we have considered the dimension 
d = 2 which is relevant for biological populations. The results of the simulations are reported 
in Figs. [D and El The uniform distribution of particles is unstable (this is similar to the Jeans 
instability in astrophysics) and the particles start to collapse and form aggregates. Note that 
the concentration can be quite large although this is not always obvious on the figures since 
the particles have fallen on each other. Density contrasts are easier to see in continuous (fluid) 
models that automatically involve a coarse-graining [T5J, [37] . However, our iV-body simulations 
also show the formation of lines and filaments in the case of models with inertia (Fig. [T]) and 
the absence (or reduction) of such lines and the generation of pointwise blow-up in the case of 
overdamped systems (Fig. [2]). 



3 Conclusion 

In this paper, we have introduced general kinetic and hydrodynamic models of chemotactic ag- 
gregation. These models can be relevant to describe the organization of social insects (swarms) 
like ants or the morphogenesis of biological populations like bacteria, amoebae, endothelial 
cells etc. Starting from a microscopic model defined in terms of iV coupled stochastic equations 
(fT3])-(fl2]) or (fl4l)-(fl5j). we have derived a generalized mean field Fokker-Planck equation ( fT7l) 
governing the evolution of the distribution function of the system in phase space. By taking 
the successive moments of this kinetic equation and closing the hierarchy by a local thermo- 
dynamic equilibrium condition, we have derived a set of hydrodynamic equations (I2"7]) - (l29]) 
involving a damping term. An interest of this approach is to connect and generalize different 
models previously introduced in the literature. In particular, the Keller-Segel model ([I])-© and 
the generalized Keller-Segel model ©-(jl]) are obtained in a limit of strong friction £ — > +00 
and the hydrodynamic model introduced by Gamba et al. ©-([7]) corresponds formally to a 
limit of low friction £ — > (in fact, the justification of this model must be given separately 
as in Appendix [D]) . We have illustrated numerically the difference between models with iner- 
tia (hyperbolic) leading to network patterns and models without inertia (parabolic) leading to 
pointwise blow-up. We have discussed the analogy between the chemotactic collapse in biology 
and the gravitational collapse (Jeans instability) in astrophysics (see also [20]). 

In our kinetic model of Sec. 12.21 the dynamical evolution of bacterial populations is described 
by generalized stochastic processes (TL4]) and nonlinear mean field Fokker-Planck equations (TT7T) 
similar to those arising in the context of generalized thermodynamics [TT] [T2] US]- As a re- 
sult, the hydrodynamic equation (125]) involves a nonlinear pressure term p(p). If we consider 
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stochastic processes with normal diffusion (IT5|) . we obtain ordinary Fokker-Planck equations 
( I66|) leading to a linear equation of state p(p) = pT e ff (isothermal) with T e ff = D/£. In the 
strong friction limit, this yields the Keller-Segel model ([1]) which is known to form singularities 
(Dirac peaks). In order to obtain the regularized model we need to modify the stochastic 
equations. This has been done here by letting the diffusion coefficient depend on the local dis- 
tribution of particles in phase space (see Eq. (JHJ)). This is a phenomelogical attempt to take 
into account complicated microscopic constraints (like excluded volume constraints, finite size 
effects, short-range interactions,...) that affect the dynamics of the particles at small scales. At 
the level of the hydrodynamic equations, these constraints are modeled by an effective pressure 
term p = p(p) which replaces the usual term p = pT e ff. An example of regularized Keller-Segel 
model where the Dirac peaks are replaced by smooth density profiles (aggregates) has been 
studied in [14] . In that case, the effective equation of state is p(p) = — &oT e ff ln(l — p/er ). 
For dilute systems where the motion of an individual cell is not impeded by the other cells, 
we have p <C cr and we recover the "isothermal" equation of state p = pT e ff. However, mod- 
ifications arise when the cells are compressed. Indeed, the equation of state departs from the 
isothermal one when the density approaches the maximum allowable density <j . In that case, 
instead of Dirac peaks, we form flat cores with density p ~ er . The nonlinear pressure term 
p(p) can also model a process of anomalous diffusion. This is the case in particular for the 
stochastic model fl32|) leading to a polytropic equation of state p = Kp 1 . This model can take 
into account effects of non-ergodicity and nonextensivity. This leads to the Tsallis statistics 
that arises when the dynamics has a fractal or multi-fractal phase space structure [2E]- The 
corresponding generalized Keller-Segel model (J4"4"|) has been studied in detail in [T71 [TSJ E5] • 

The next step is a detailed study of the models presented in this paper. This study is of 
interest not only in mathematical biology but also, at a more general level, in statistical me- 
chanics. What we are studying essentially is the Dynamics and Thermodynamics of Brownian 
Particles with Long-Range Interactions [2B]. We have furthermore introduced a general class 
of stochastic equations with a random force depending on the distribution function, forcing 
the nonlinearity in the Fokker-Planck equation. This is associated with a notion of general- 
ized thermodynamics. Therefore, our model involves both long-range forces pQ and generalized 
thermodynamics (related to nonlinear Fokker-Planck equations) [19], which are two domains 
actively studied at the moment in statistical mechanics. Our approach shows that these topics 
can have applications in biology. In the overdamped limit, the dynamical equations have the 
form of drift-diffusion equations corresponding to the Keller-Segel model or generalizations of 
this model. These equations have been extensively studied in the mathematical (see the re- 
view of Horstmann [5]) and physical (see Chavanis & Sire [21] and references therein) literature. 
The more complicated study of kinetic equations (IT7|) - (IT5|) or hydrodynamic equations (I2"7j) -(l2"9"l) 
should be considered in future works. The connection with generalized forms of Cahn-Hilliard 
equations when the potential of interaction is short-ranged (as noticed in [T^l E3J) should also 
be developed (see Appendix [F]). 

A Many-body theory of chemotactic aggregation 

In this Appendix, we first develop an exact many-body theory of Brownian particles in interac- 
tion. Then, we show how a mean field approximation can be implemented in the problem. For 
simplicity, we assume that each particle has a normal Brownian motion, i.e. the diffusion coef- 
ficient D is constant. Our approach follows the steps of Newman & Grima [22]. However, we 
take into account the inertia of the particles while Newman & Grima consider an overdamped 
limit. 
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Basically, the dynamical evolution of iV Brownian particles in interaction is described by 
stochastic Langevin equations of the form (fT3|) . In the biological context, the field c represents 
the concentration of the chemical produced by the cells and it satisfies an equation of the form 
(|T2|) . We define the single-cell probability distribution in phase space by 



P«(r,v,t) = (5{r-r a {t))S{v-v a {t))), 



(54) 



where the brackets denote an average over the noise. Similarly the two-cell probability distri- 
bution in phase space is 

P^(r,v,t;r',v',0 = (5(r - r a {t))5{v - v a {t))8{r' - r p (t'))8(v' - v^i')))- (55) 

Integrating Eq. (1121) . the concentration field can be expressed in terms of the cell paths as 

c(r, t) = h [ dr' [ dt'Giv -r',t- t') V 8{r' - r a (t')), 
J Jo 

where the Green function for the chemical diffusion equation is 



(56) 



G{r,t) = (4itD c t)- d/2 exp 



^D c t 



- kt 



(57) 



Taking the time derivative of P a , we get 



8P 8 8 

-^ = -JT- " r«(*))*(v - v a (t))> - j- ■ {v a (t)S(r - r a (f))*(v - v a (t))>. (58) 

8t or C7v 

Inserting the equations of motion ffl3l) in Eq. (1581) . we obtain 
BP 8 8 

-jr- ■ (V a c 6(r - r a (t))<5(v - v a (*))) - tt" ■ (v^R^t) 5(r - r a (t))5(v - v a (t))>. (59) 
av av 

The first two terms are straightforward to evaluate. The fourth term is the standard term that 
appears in deriving the Fokker-Planck equation for a pure random walk; it leads to a term 
proportional to the Laplacian of P a in velocity space. The third term can be evaluated by 
inserting the formal solution f|56l) in Eq. (1591) . This leads to the exact equation 



BP a 

8t 



+ v-^ + /i^-- [dr 1 [ dt'VG(r-r',t-t')YP a 
dr 8v J J Q *-f 



p{r,v,t; r',t') 



8_ 
8v 



8P 

D^ + ZP a v 
av 



where the statistical correlations are encapsulated in the two-cell distribution 
P^(r,v,t;r',t') = (8(r - r„(f))<f(v - v a (t))S(r' - r p {t'))). 



(60) 



(61) 



The mean field approximation is implemented by assuming that the two-cell distribution fac- 
torizes 



P a> p{r, v, t; r', t') = P a {r, v, t)P p {r', t'). 



(62) 
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In that case, the preceding equation can be rewritten 
dp no a r Ft 



dt 



v-^ + ^-P Q (r,v,t)V fdx'f d£G(T-T?,t-H)Y,Pf>VJ 

j3 



d_ 
dv 

We define the distribution function and the density by 



N 



a=l 



BP 
ov 



(63) 



/(r,v,t) = ^P Q (r,v,t), 



(64) 



N 

p(r,t)=^P a (r,t) = f(r,v,t)dv. 

a=l J 



(65) 



Summing Eq. fl63|) over the cell index, we obtain a kinetic equation of the form 



dt dr <9v dv 



<7V 



where 



c(r,t) = /i y rfr'^ <tfG(T-T , ,t-t')p(i',t). 



This smooth field satisfies the partial differential equation 

<9r: 



dt 



—kc + D c Ac + hp. 



(66) 



(67) 



(6* 



Note that Eq. (166]) can be viewed as a mean field Kramers equation [26]. For £) = £ = 0, it 
becomes equivalent to the Vlasov equation. 

If we consider a limit of strong frictions (or large times t 3> £ -1 ) in the stochastic equations 
(1T3"]) . we obtain the overdamped model (ITT]) . These equations, coupled to Eq. (|T2"|) for the 
chemical concentration field have been considered in [22] • The above procedure leads, in the 
mean field approximation, to a kinetic equation of the form ([1]) coupled to Eq. fl68l) . This returns 
the ordinary Keller-Segel model. Equation can be viewed as a mean field Smoluchowski 
equation [26]. It can also be obtained from the mean field Kramers equation Eq. (166]) by 
considering the strong friction limit and using a Chapman- Enskog expansion [53] • Therefore, 
our model (fT3l can be considered as an extension of the overdamped model of [23] when inertial 
effects are taken into account. Finally, inertial and overdamped stochastic models of the form 
(TLB"]) and (TTTT) where the particles interact via a binary potential (instead of a field equation (fl2l) 
depending on the history of the system) have been studied in [26J. A hierarchy of equations 
has been derived for the reduced probability distributions and Eqs. (1561) and (JT]) are obtained 
in a mean field approximation valid in a proper thermodynamic limit with iV — > +oo. 

Let us finally consider the deterministic case without diffusion D = T e ff = 0. We introduce 
the exact distribution function and the exact density 



N 



f ex (r, v,t) = £ 6{r - r a (t))«J(v - v Q (t)), 



(69) 



a=l 
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N 

Pex {r, t) = J2 S(t - r a {t)) = / f ex (r, v, t)dv. (70) 

a=l J 

Then, from the equations of motion (TT3T) . repeating the steps (T58l - (T60l) without the brackets, 
we obtain the exact equations 

^T + v-^r + Vc,, — = -.(^/ e ,v), (71) 

dc 

= -kc ex + D c Ac ex + hp ex . (72) 

These equations bear exactly the same information as the iV-body system (fTBI) with D — 0. 
For £ = 0, Eq. (17T|) becomes equivalent to the Klimontovich equation in plasma physics [38J. 
On the other hand, in the strong friction limit, if we start from the iV-body system ( TTTT) with 
= T e ff = 0, we obtain the exact equations 

dpex 



dt 



), (73) 
-kc ex + D c Ac ex + hp ex . (74) 



dt 

These equations bear exactly the same information as the iV-body system (TTTT) with .D* = 



B The limit of strong friction 



The generalized Smoluchowski equation ( |31j) can be derived from the generalized Kramers 
equation (TTTT) in the strong friction limit £ — > +oo. Using the effective Einstein relation £ = 
the generalized Kramers equation (TTTj) can be rewritten 



^ + v 21 + V5( c ) — - — 



fC"(f)^-+i3fv 



(75) 



In the strong friction limit £ = D/3 — ► +oo, assuming f3 of order unity, the term in bracket in 
Eq. (1751) must vanish (since -D — ► +oo) so that the distribution function satisfies to leading 
order 



y + A(r,t) 



0(£ 



(76) 



C7'(/) = -/3 

The function A(r, t) is related to the spatial density p(r,t) through the relation 

p(r,t)= j fdv = p[\(v,t)\. (77) 

Using Eq. (1761) . we find that u = 0(£ _1 ) and = p5ij + 0(£ -1 ) where p(r, t) is the local 
pressure 

p(p,t) = i | /w a dv = p[A(r,f)], (78) 

determined from Eqs. (T76|) and (TTTj) . As in Sec. 12.31 the fluid is barotropic, i.e. p(r, t) = 
p[p(r, t)], where the equation of state p(p) is obtained by eliminating A(r,t) between Eqs. ( 1771) 
and (1781) . The equation of state is entirely specified by the function C(f). To first order in 
the momentum equation (12~T1) implies that 

pu = -I(Vp - pS'(c)Vc) + 0(C 2 )- (79) 
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Inserting the relation (1791) in the continuity equation (T2"Uj) . we get the generalized Smoluchowski 
equation (l3Tj) . The generalized Smoluchowski equation, as well as the first order correction 
to the distribution function /(r, v, £), can also be obtained from a formal Chapman-Enskog 
expansion [33]. Note finally that, instead of the generalized Fokker-Planck equation (175|) . we 
could consider the generalized isotropic BGK equation 



df df df f - f 

- + v.- + V5(c).- = -^, (80) 

where f (r,v,t) is defined by C'(fo) = —f3(v 2 /2 + A(r, t)) where A(r, t) is determined by 
J fodv = p(r,t). The calculations are similar to those described previously with 1/r play- 
ing the role of £. In particular, the generalized Smoluchowski equation (13T|) is obtained at order 
0(t) for r — > 0. A Chapman-Enskog expansion is given in Appendix A of [33J. 

C The limit of large difFusivity 

Let us consider the Keller-Segel model 

^ = J D,Ap- X V-(pVc), (81) 
dc 

— = -kc + hp + D c Ac, (82) 
with Neumann boundary conditions 

Vp • n = Vc • n = 0, (83) 

where n is a unit vector normal to the boundary of the domain. Following Jager & Luckhaus 
[39J, we set h = \D C and assume that D c — > +oo and A ~ 1. Introducing the average density 
p = (1/V) / pair (where V is the volume of the box), we find from Eq. (18T|) that p(t) = p . 
Then, from Eq. (1521) . we get 

ft + ^ = APo- (84) 



We note, parenthetically, that this equation has the exact solution, 

c(t) = ^p (l-e- fct )+c e- fci , (85) 

so that the average concentration of the chemical relaxes to c(+oo) = (h/k)p on a timescale 
k~ l . If we consider c = c — c, we find from Eqs. (l82l) and ( |84|) that 

1 ( dc \ 

kc) = Ad + X(p-p ). (86) 



In the limit D c — > +oo, we obtain the reduced Keller-Segel model 

^ = J D,Ap- x V-(pVc), (87) 

Ac = -A(p-p ). (88) 

Note that this system is mathematically well-posed with the Neumann boundary conditions 
(183|) . contrary to the case where Eq. (1881) is replaced by the Poisson equation 

Ac = -Ap. (89) 
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Indeed, if we integrate Eq. (|89|) over the domain and use the divergence theorem, we obtain a 
contradiction as 

r Vc • dS = [ Ac dr = -AM ^ 0. (90) 



However, when the density blows up so that p ^> p , it is justified to consider the model 
(|87p - (l89p as an approximation. It is only in that case (large diffusivity of the chemical and 
high concentration of the particles) that the Keller-Segel model for the chemotaxis becomes 
equivalent to the Smoluchowski-Poisson system for self-gravitating Brownian particles. On the 
other hand, we note that a homogeneous distribution is an exact stationary solution of the 
Keller-Segel model (|81|) - (|82p with kc = hp. It is also an exact stationary solution of the 
reduced Keller-Segel model (j8"7j) - (j8"8"l) with p = ~p . By contrast, a homogeneous distribution is 
not a stationary solution of Eqs. (1571) and (I89I) when Eq. (1881) is replaced by a Poisson equation 
( |89|) as in astrophysics. In astrophysics, we have to advocate the Jeans swindle when we analyse 
the linear dynamical stability of an infinite and homogeneous self-gravitating system [3D]. By 
contrast, there is no "Jeans swindle" in the chemotactic problem if we properly use Eqs. ( 1821) 
or (EED instead of Eq. (|89l) [20]. 

Setting k = k^D c and h = XD C , we can also consider the limit D c — > +oo with A, k ~ 1. 
Equation (|82|) can be rewritten 

1 dc 

— -^ = -k 2 c + Xp + Ac, (91) 

which reduces, for D c — > +oo, to 

Ac-klc=-\p. (92) 



D Kinetic derivation of the hyperbolic model 

The kinetic theory presented in Sec. [2] is well-suited to weakly inertial systems for which the 
damping term — £u is relatively strong. For £ — > +oo, we rigorously obtain the generalized 
Smoluchowski equation ( }3~TT) . For large values of £, the damped Euler equations (J27j)-(l29l) may 
provide a relatively good description of the dynamics (but we again stress that they are not 
rigorously justified). Alternatively, the model ©-© proposed by Gamba et al. [15] corresponds 
to highly inertial systems. Formally, it can be obtained from Eqs. (|27|) - (|29|) by taking — £u = 0. 
However, we have indicated in Sec. 12.31 the shortcoming of this procedure when we start from 
a stochastic equation of the form ( fl3|) . Indeed, the L.T.E. condition ( 1241) is not justified when 
£ — * 0. Here, we present an alternative kinetic model leading rigorously to the model ©-(CO) 
proposed by Gamba et al. [15J. We assume that the particles obey a stochastic equation of the 
form 

dv 



(Jf -C(v - u(r, *)) + VS(c) + v / 2 J D(/)R(t). (93) 

The first term in the r.h.s. is a friction force relative to the mean velocity u(r, t). Its physical 
effect is therefore completely different from the friction force in Eq. ( fl4l) . The associated 
nonlinear mean field Fokker-Planck equation is 



?L + v M + vs(c).?t = — .tD 
dt dr dv dv 



fC"(f)^ + Pf(v-u(v,t)) 



(94) 



where (3 = C/D. Taking the hydrodynamical moments of this equation, as in Sec. 12. 3[ we 
obtain the continuity equation (1201) and 

d d dP dc 

M {fMi) + d^ {pUiU ^ = ~^; + ^' (C W (95) 
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Contrary to Eq. (|21|) . this equation has no damping term and it conserves the total impulse. 
Considering the limit ( = Df3 — > +00 with (3 fixed, we find from Eq. (1941 that the distribution 
function satisfies 



C'tf) = -(3 



w 2 



A(r,t) 



OiC 1 ). (96) 



Therefore, in this approach, the L.T.E. (1241) is exact to leading order in ( — > +00. This implies 
that Pij = p5ij + 0(C _1 ) where the pressure p(r,t) satisfies a barotropic equation of state p(p) 
determined by the function C(f) as in Sec. 12.31 Substituting this result in Eq. fl95l) . we get the 
hydrodynamic model ©-([7]) without friction term proposed by Gamba et al. [15]. Note finally 
that, instead of the generalized Fokker-Planck equation fl94l we could consider the generalized 
BGK equation 

df df df f - f l.t.e. , nvA 
w + v.- + V5(c).- = , (97) 

where /l.t.b.(i", v, t) is defined by Eq. ( l24l) . The calculations are similar to those developed 
previously with 1/r playing the role of (. In particular, the model ©-(Ej) is obtained at order 
0(t) for r — > 0. Note, finally, that other kinetic theories for chemosensitive movement have 
been proposed in [37J. They lead to hydrodynamic equations of the form (J5j)-([TJ) , without 
friction force. 



E Generalized free energies and Lyapunov functionals 

In this Appendix, we show that the different equations introduced in this paper are associated 
with Lyapunov functionals that can be interpreted as effective "generalized free energies" . 



E.l Kinetic model 

Let us first consider the kinetic model (fT71)-(|T5l) with constant coefficients 





dt dr <9v <9v 



dc 
~dt 



DfC"(f)^ + £fv 



-kc + hp + D c Ac. 
Recalling that T e ff = D/£, we introduce the functional 

F[f] = J fY drdv+ ^h J i D ^ c ) 2 + kc2 ] dr ~ J pcdr + T eff J C{f)dvdw. 



After simple algebra, one can show that 



J jj U>fC"{f)^ + £/v) drdv - i J (D c Ac -kc + hpfdv < 0. 



(98) 
(99) 

(100) 
(101) 



Therefore, F(t) decreases monotonically with time and plays the role of a Lyapunov functional. 
It can also be interpreted as a generalized free energy [ 1 1 J . This is particularly clear if we 
consider the field equation 

D c Ac -kc= -hp, (102) 
instead of Eq. (j99"]) [see Appendix [U]. In that case, the functional (HOOP reduces to 

.2 



/ — drdv / pcdr + T e ff 



C{f)dvdv. 



(103) 
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This can be written F = E — T e ffS where E = J f^drdv — | J pcdr = K + W represents the 
energy (kinetic + potential) and S = — J C(f)drdv a generalized entropy. In that case, we 
have 

F = - j jj (DfC"(f)^ + £/v) drdw < 0. (104) 

This inequality can be viewed as an appropriate if -theorem in the canonical ensemble [TTJ [26] , 
where the effective temperature T e ff is fixed instead of the energy In the case of normal 
diffusion C(f) = / In/, the functional (11031) coincides with the Boltzmann free energy 

F = J f-drdw - l -J pcdv + T eff J fin fdrdv. (105) 

We also note, parenthetically, that in the absence of diffusion (D = T e ff = 0), the free energy 
reduces to the energy. Therefore, the proper if -theorem becomes E < 0. More precisely, we 
have E = — £ j fv 2 dvd\ = —2^K < 0, where K is the kinetic energy. 

The steady state of the system (I98l) -(j99l) must satisfy F = 0. According to Eq. (11011) . this 
implies that the current in Eq. (1981) vanishes 

C"(/)|£ + 0v = O. (106) 

Then, the condition df/dt = implies that the advective term must also vanish 

v-^ + Vc-^ = 0. (107) 
or av 

Therefore, the advective term and the current vanish independently. Equation (11061) can be 
integrated on the velocity to yield 

C"(/) = -/3y-A(r), (108) 

where A(r) is an arbitrary function of the position. Differentiating this expression with respect 
to r and v, we obtain 

C"U)% = -VA(r), C"(f)?f = -(3v. (109) 

or (TV 

Substituting these relations in Eq. fllQTj) . we get 

(VA + pVc) ■ v = 0, (110) 
which must be valid for all v. This yields V 'A + /3Vc = 0, so that 

A(r) = -/3c(r) + a, (111) 

where a is a constant of integration. Substituting this result in Eq. ( 11081) . we obtain the steady 
state (|T8l) . Therefore, at equilibrium, the distribution function depends only on the energy: 
/ = /(e) with e = v 2 /2 — c(r). This first implies that the local velocity u(r) = 0. The density 
and the pressure can then be written p = J f(e)dv and p — \ J f(e)v 2 dv. Since 

Vp = -ivc j f'(e)v 2 dv = -^Vc J ^- ■ vdv = Vc J fdv = pVc, (112) 
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we find that the condition f — f(e) implies the condition of hydrostatic equilibrium 



Vp-pVc = 0. (113) 

We can also justify the Local Thermodynamical Equilibrium (LTE) assumption made in 
Sec. 12.31 from the free energy (llOOp . The idea is to close the hierarchy of hydrodynamic 
equations by calculating the pressure tensor (!22j) with the distribution /lte that minimizes the 
free energy (1 1001) at fixed p(r, t) and u(r, t). If we prescribe the density p(r, t), the concentration 
of the chemical c(r, t) is automatically fixed by Eq. (1991) . This implies that the second and 
third integrals in Eq. (llOOp are fixed. Therefore, minimizing F at fixed p(r, t) and u(r, t) is 
equivalent to minimizing 

fldvdw + T eff J C{f) drdv, (114) 

at fixed p(r, t) and u(r, t). Introducing Lagrange multipliers and writing the variational problem 
in the form 

5F - J a(r, t)5f drdw -J b(r,t)5f-v drdw = 0, (115) 

we obtain 



v 2 



- + T eff C'(f) - a(v,t) - b(r,t) ■ v = 0. (116) 

Relating the Lagrange multipliers a(r,t) and b(r, t) to the constraints p(r,t) and u(r, t), we 
obtain the distribution ([24]). Since 5 2 F = (l/2)T e// / C"(f)(5f) 2 drdv > 0, and C is convex, 
the distribution (j24"|) is a minimum of F at fixed p(r,t) and u(r, t). 

E.2 Fluid model 

Let us consider the damped barotropic Euler equations with constant coefficients 

^ + V-(pu) = 0, (117) 

— + (u-V)u= — Vp + Vc-^u, (118) 
ot p 

dc 

— = -kc + hp + D c Ac. (119) 
at 

We introduce the functional 

F[p, u] = J p J" ^p-dp'dv [D c (Vc) 2 + kc 2 } dr- J pcdr + pu 2 dr. (120) 

This functional can be obtained from Eq. (llOOj) by using the LTE condition fl24l) to express 
F[f] as a functional of p and u (the calculations are similar to those detailed in [T3| 33J). After 
simple algebra, it is easy to show that 

F = -£ I pu 2 dr - - I (D c Ac - kc + hpfdr < 0. (121) 
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Therefore, F(t) decreases monotonically with time and plays the role of a Lyapunov functional. 
It can also be interpreted as a generalized free energy. In the case where Eq. fl 1 1 9 [) is replaced 
by Eq. (jID2l> . we have 

F[p, u] = J pJ P ^dp'dr - l -j pc dv+ l -j pu 2 dr. (122) 

and 

F = -f / pu 2 dr < 0. (123) 



Note that for an isothermal equation of state p = pT e ff, corresponding to normal diffusion, the 
functional (11221) returns the Boltzmann free energy 

F[p,u] =T e ff J plnpdr-^ J pcdr + ^ J pu 2 dr. (124) 

The steady states of Eqs. (1 11 7ft - (1 11 91) must satisfy F = yielding, according to Eq. (I12ip . 
u = 0. Then, using Eqs. (I117p - (I119I) . we obtain the condition of hydrostatic equilibrium flll3j) . 

E.3 Overdamped model 

Let us consider the generalized Smoluchowski equation with constant coefficients 

(125) 



dt 



- (Vp - pVc) 



3c 

— = -kc + hp + D c Ac. (126) 



We introduce the functional 



F[p] = J pf ^dp'dv + i- J [D c (Vcf + kc 2 } dr- J pcdv. (127) 

This functional can be obtained from Eq. (11001) by using Eq. f!76|) to express F[f] as a functional 
of p in the strong friction limit [TBI [33] . It can also be obtained from Eq. (I120p by using the 
fact that u = 0(^ _1 ). After simple algebra, it is easy to show that 

F — — f —(Vp- pVcfdr - - I (D c Ac -kc + hpfdv < 0. (128) 
J SP h J 

Therefore, F(t) decreases monotonically with time and plays the role of a Lyapunov functional. 
It can also be interpreted as a generalized free energy. In the case where Eq. f)126p is replaced 
by Eq. ffTU2D . we have 

F\p\ = J pf yrdP^ ~\J P^r, (129) 

and 

F — — —(Vp- pVcfdv < 0. (130) 

J ip 

Note that for an isothermal equation of state, p = pT e ff, we obtain the Boltzmann free energy 
in configuration space 

F[p]=T eff J plnpdr-^ J pcdv. (131) 

The steady states of Eqs. f ll25l) - fll26l) must satisfy F = yielding, according to Eq. (11281) . the 
condition of hydrostatic equilibrium (I113p . 
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E.4 The primitive Keller-Segel model 

The primitive Keller-Segel model has the form [I]: 

^ = V • (D 2 Vp) - V • (DiVc), (132) 
Oc 

— = -k(c)c + h(c)p + D c Ac, (133) 

where Z?i = Di(p,c) and L> 2 = D 2 (p,c) can both depend on the concentration of cells and of 
the chemical. Let us consider a simplification where D 2 = Dh(p), D\ = X9{p)i M c ) = & an d 
/i(c) = h. In that case, we obtain 

^ = V • [Dh(p) Vp - Vc] , (134) 
= _jfec + /i P + £) c Ac. (135) 

at 

Equation (11341) can be viewed as a nonlinear Fokker-Planck equation of the form considered in 
[TTj . It can be obtained from the stochastic process 

dr 



- = X (p)Vc+y/2D(p)K(t) } (136) 

where R(i) is a white noise, x{p) = X9{p)/p an d -^(p) = {D/p) J P h(p')dp'. When g(p) = p and 
h(p) = 1, we recover the ordinary Keller-Segel model (pQ) with constant diffusion coefficient D 
and constant mobility x- When g(p) = p and h(p) is arbitrary, Eq. (11341) describes a situation 
where the mobility is constant but the diffusion coefficient can depend on the density. This 
can account for anomalous diffusion and non-ergodic effects. This is the case in Eqs. (1441) and 
( |4"T|) corresponding to g(p) = p and h(p) = (x/ D)K^p y ~ l [17J or, more generally, in Eq. (j3J) 
corresponding to g(p) = p and h(p) = (x/D)p'(p) [TT]. When h(p) = 1 and g(p) is nonlinear, 
Eq. (11341) describes a situation where the diffusion coefficient is constant but the mobility 
(or chemotactic sensitivity) depends on the density. This can account for excluded volume 
effects or steric hindrance when the density is high. For example, the case where h(p) = 1 
and g(p) = p(l — p/a ) has been treated in [HI [TJ]. In that case, the mobility tends to zero 
when the density approaches its maximum value <Jq. More generally, the drift-diffusion equation 
(11341) describes a situation where the diffusion coefficient D(p) and the mobility x(p) can both 
depend on the density. It can be derived from a master equation by allowing the transition 
probabilities from one site to the other to depend on the occupancy number [30j IT4"] . 
If we introduce the functional 

F[p] = - J C(p) dv + i- J [D c (Vc) a + kc 2 ] dv - J pcdr, (137) 

where C(p) satisfies C"(p) = h(p)/g(p), it is easy to show after simple algebra that 

F — — / —^—{Dh(p)Vp - xgi.p)^c) 2 dv - - I (D c Ac - kc + hpfdv < 0. (138) 
J X9{P) h J 

Therefore, F(t) decreases monotonically with time and plays the role of a Lyapunov functional. 
In the case where Eq. f l 1 3 5 [) is replaced by Eq. (11021) . we have 

F[ P } = ^jc(p)dr- 1 - jpcdr, (139) 
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and 



F = - / -— -(D%)Vp - X g{p)Vcydr < 0. (140) 
J X9{P) 

This can be written F = E—T e ffS where E = —(1/2) J pcdv is the potential energy, T e ff = 1/13 
is an effective inverse temperature satisfying a generalized Einstein relation (3 = xjF> and 
S = — J C(p)dr is a generalized entropy. Therefore, F can be interpreted as a generalized free 
energy in an effective thermodynamical formalism. We can thus write Eq. (11341) in different 
forms as shown explicitly in [TTj [T4]. This strengthen the analogy with generalized Fokker- 
Planck equations. The steady states of Eqs. fl 1 341) - (1 1 3 5 [) must satisfy F = 0. According to Eq. 
(11381) . this implies that the current must vanish yielding the relation 

C'( Peq )=(3c-a, (141) 

where a is a constant of integration. Since C is convex, the above relation can be reversed. 
Then, we find that p e g( r ) = F(—(3c + a) = p eq (c) where F(x) = (C") _1 (— x). Since p' eq (c) = 
(3/C"(p eq ) and (3 > 0, we find that p eq (c) is a monotonically increasing function of the con- 
centration. Equation (11411) can also be obtained by extremizing the free energy (11371) at fixed 
mass. Furthermore, it is shown in [TT] that linearly dynamically stable solutions of (1134^ - (1135H 
correspond to minima of F at fixed mass (these general results also apply to the other model 
equations discussed previously). Finally, we note that Eq. f |134j) can be written in the form 



dt 



(142) 



where S/Sp denotes the functional derivative [12J. 



F Generalized Cahn-Hilliard equations 

In this Appendix, we show that, in the limit of short-range interactions, the kinetic equations 
presented in this paper reduce to generalized forms of the Cahn-Hilliard equation describing 
phase ordering kinetics [42]. The connection to the Cahn-Hilliard equation was previously 
mentioned in [T2| 133]. 

Let us first consider the situation where the equation for the chemical is given by (see 
Appendix 

Ac- klc= -Xp, (143) 

In the limit ko — > +oo, we can neglect the Laplacian and obtain in first approximation c ~ 
Ap/fep. Then, substituting this relation in the Laplacian we get the next order correction 

c^hAp + k 2 p). (144) 

More generally, suppose that the concentration of the chemical is determined by a relation of 
the form 

c(r,t) = J u{\r-r'\)p(r',t)dr', (145) 

where u is a binary potential of interaction. Let us now assume that u(|r — r'|) is a short-range 
potential of interaction. Then, setting q = r' — r and writing 

c(r,t) = / u(g)p(r + q,t)rfq, (146) 
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we can Taylor expand p(r + q, t) to second order in q so that 



p(r + q) t) = p(r,t) + ^^ + i^^ 



Substituting this expansion in Eq. (11461) . we obtain 



dxidx 



-qiQj- 



(147) 



c(r,t) = ap(r,t) + -Ap(r,t), 



with 



Sri 



+ 00 



u(q)q d 1 dq 



d 



s d 



+00 



u{q)q d+1 dq. 



:i48) 



(149) 



Note that I = (b/a) 1 ^ 2 has the dimension of a length. In the limit of short-range interactions, 
we can replace c(r, t) by Eq. (11481) in all the dynamical models introduced in this paper. Of 
particular interest is the case of the generalized drift-diffusion equation (1 134ft . Substituting Eq. 
(11481) in Eq. (1 134p and introducing the effective potential 



V( P ) 



a o 2£> 

._„ 2 + _ C(/ ,) + Vi, 



we get 



^ = -AV-[g(p)V(Ap-V'(p))} 



(150) 



(151) 



with A = ybjl. Morphologically, this equation is similar to the Cahn-Hilliard equation |42j. 
One difference, however, is that g(p) can depend on the density (for constant mobility g(p) = p) 
while this quantity is constant in the ordinary Cahn-Hilliard equation. With this (important) 
difference in mind, the drift-diffusion equation (11341) can be seen as a generalization of the 
Cahn-Hilliard equation to the case of long-range potentials of interaction. For the particular 
case D = T e ff = 0, we get 



dp 
dt 



-AV ■ 



(/(/>) V ( Ap+^-p 



fl52) 



Note, finally, that Eq. (1151R can be written in the form 



where 



F[p] 



dr. 



(153) 



(154) 



This expression of the free energy can be obtained from Eq. fl 1 3 9 j) . by using Eq. (j!48p . For 
comparison, the ordinary Cahn-Hilliard equation for model B (conserved dynamics) is [42J: 



dp 
di 



A 



6£ 



(155) 



In the Cahn-Hilliard problem, the potential has a double-well shape leading to a phase sep- 
aration while, in the present case, the potential can take other forms, as in Eq. (11521) for 
example. 
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